Replacing standard galaxy profiles with 
mixtures of Gaussians 
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ABSTRACT 

Exponential, de Vaucouleurs, and Sersic profiles are simple and successful 
models for fitting two-dimensional images of galaxies. One numerical issue en- 
countered in this kind of fitting is the pixel rendering and convolution (or corre- 
lation) of the models with the telescope point-spread function (PSF); these oper- 
ations are slow, and easy to get slightly wrong at small radii. Here we exploit the 
realization that these models can be approximated to arbitrary accuracy with 
a mixture (linear superposition) of two-dimensional Gaussians (MoGs). MoGs 
are fast to render and fast to affine-transform. Most importantly, if you have a 
MoG model for the pixel-convolved PSF, the PSF-convolved, affine-transformed 
galaxy models are themselves MoGs and therefore very fast to compute, inte- 
grate, and render precisely. We present worked examples that can be directly 
used in image fitting; we are using them ourselves. The MoG profiles we provide 
can be swapped in to replace the standard models in any image-fitting code; they 
sped up model fitting in our projects by an order of magnitude; they ought to 
make any code faster at essentially no cost in precision. 

Gaussians are remarkable distribution functions. They have the incredible properties 
that — in any number of dimensions — the convolution (or correlation) of one multivariate 
Gaussian with another is itself a multivariate Gaussian, and any product of multivariate 
Gaussians is itself a multivariate Gaussian but with a different normalization. Furthermore, 
the means and variance tensors of the Gaussians output by these operations are related 
simply to the means and variance tensors of the inputs. Add to these wonders the fact that 
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Gaussians form a complete basis for representing (smooth) probability distribution functions 
and it becomes remarkable that we don't do everything we do in terms of Gaussians. 

To elaborate, a mixture of multivariate Gaussians — a linear superposition — can be used 
to represent any reasonable distribution in any number of dimensions to any reasonable pre- 
cision. Convolution (or correlation) by any other distribution that has also been represented 
by a mixture of multivariate Gaussians creates a new mixture of Gaussians with simply 
adjusted amplitudes, means, and variance tensors. The ubiquity of convolution operations 
in astronomy suggests the widespread adoption of mixture-of- Gaussian (MoG) modeling. 
We have pioneered this in the area of distribution modeling in large numbers of dimensions 
(Bovy etal. 2011a|[b 2012), where convolution occurs because the true or noise-free distri- 
bution is convolved with the noise before being observed. Here we are going to capitalize on 
the convolution properties of MoGs in modeling galaxy morphologies in imaging data, where 
the true or high-angular-resolution intensity field is convolved with the point-spread func- 
tion (PSF) before being observed. This is always easier with Gaussians than cuspy profiles, 
but is particularly useful in data sets in which the PSF itself has also been modeled as a 
MoG, which is not uncommon (for example, the Sloan Digital Sky Survey imaging pipelines 
described by Lupton e^a/.|200T make approximate MoG PSF models for every imaging field). 



We are not the first in this space: Deconvolution and modeling of galaxy images with 



MoGs has been done very successfully before (for example, Bendinelli 1991 Emsellem etal 



1994 Bendinelli &; Parmeggiani 1995 |Cappellari |2002[ they use the name "multi-Gaussian 
expansion" or MGE). However, the idea behind those projects was to use the MoGs to provide 
a very free form for the isophotes or morphologies of resolved galaxies or other complex 
scenes. Here our goals are very limited: We want to improve the performance of standard 
galaxy image model fitting by expressing the standard galaxy models — the exponential and 
de Vaucouleurs profiles — as rigid MoGs. 

Whenever an investigator is fitting PSF-convolved exponential or de Vaucouleurs pro- 
files, the models presented here will improve code performance. That doesn't mean that 
doing such fitting is a Good Idea. These profiles are effective models for galaxies at low 
signal-to-noise (as they are used in the SDSS) and they are useful templates for perform- 



ing consistent photometry across morphologically diverse objects (as we use them in Bundy 



etal. 2012). More radically, these profiles are sometimes used to distinguish disk and bulge 



light (as in Simard et al. 2011 and references cited therein); these uses are prone to over- 
interpretation: There is no theoretical argument and only weak observational arguments that 
the rotation-supported parts of galaxies are always exponential and that the kinematically 
hot components are always more de-Vaucouleurs-like. 



The performance advantages we will obtain from using MoG approximations are not 
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simply that convolution itself is trivial. The standard unconvolved galaxy models, espe- 



cially the de Vaucouleurs and Sersic models (de Vaucouleurs 1948, Sersic 1963), are very 



ill-behaved near the galaxy center. Rendering these profiles precisely near the center can be 
very challenging numerically. In addition, most of the rendering time is spent at very small 
radii, where PSF convolution is going to erase all structure anyway. That is, a MoG descrip- 
tion of the de Vaucouleurs profile saves time both in rendering and in PSF convolution; it 
produces profiles that are approximate but very high performance when all real uses of the 
profiles are PSF-convolved, as they usually are. 

On the subject of PSF-convolution, it is important in any image-modeling situation to 
think of the PSF as the pixel- convolved point-spread function. Under this choice, synthesis of 
a pixelized image involves only convolution with the PSF and evaluation at the pixel centers. 
That is, image synthesis (returned to at the end of this Note) consists of PSF convolution — 
an arithmetic operation on the two MoGs — followed by evaluation of the resulting MoG at 
the pixel centers). 

We need good performance in two-dimensional image synthesis (for fitting) because, 
with The Tractor (Lang etai, forthcoming), we are building a comprehensive model of all 
the imaging data we have; this will have models of many millions of galaxies in imaging 
that contains on the order of 10^^ pixels. We are basing our models on the Sloan Digital 
Sky Survey (SDSS) Catalog galaxy models, which only include exponential, de Vaucouleurs, 
and composite (mixture of the two) radial profiles. In detail, in fact, the SDSS Catalog 
models are small modifications of these; details below. So our goal here is to provide re- 
placements for these models in order to improve the performance of galaxy image modeling 
and analysis software. Some of the models we present have been used previously in tools 



for precise photometry (Bundy etal. 2012), and are being used in the Large Synoptic Sur- 



vey Telescope (LSST) prototype galaxy photometry pipeline (Shaw 2012), in which these 



models are convolved with a shapelet representation of the PSF; because shapelets are sim- 
ply polynomial perturbations on a Gaussian, they can also be convolved analytically with 
mixture-of-Gaussian galaxy models ( Bosch|[2010 ). These models could also be used to speed 



other image-fitting systems, like the very successful GALFIT (Peng etal. 2002). An al- 



ternative to making MoG profiles is to find exact analytic expressions for certain kinds of 
convolutions. There are some analytic results for PSF-convolved Sersic models but they in- 



volve some non-trivial series and special-function expressions (Trujillo etal. 2001); we won't 



come back to these again but they might be very useful in many situations. 

Because we are thinking about two-dimensional imaging, we use here two-dimensional 
Gaussian or Normal distributions, which look like 
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where x and m are two-dimensional vectors (usually in the focal plane or on the sky or 
something like that), and V is a symmetric 2x2 variance tensor or matrix, and implicitly 
the vectors are column vectors. A mixture of Gaussians is a linear superposition of Gaussians. 
Any positive two-dimensional function with finite support and finite total integral — including 
as a special case any two-dimensional probability distribution function — can be represented 
as a MoG to arbitrary accuracy; that is 

K 

p{x) ^ ^akN{x\nik,Vk) (2) 

k=l 
K 

1 = ^Ofe , (3) 

k=l 

where p{x) is any probability distribution of a two-dimensional quantity x, the ~ symbol 
implies approximation, K is the number of Gaussians used in the MoG, and the K Gaussians 
have amplitudes a^, means m^, and variance tensors Vk- The sum-to-one condition ensures 
that the probability distribution approximation is properly normalized. 

We wish to make an approximation to the two-dimensional circular exponential (exp) 
profile g<^^P(-) 

Q-P(C) = exp(-a-P[|C|-l]) (4) 
a^^P = 1.67834699 , (5) 

where £ is a dimensionless focal-plane position, and is a dimensionless inverse length set 
to ensure that the profile has unit half-light radius. The position ^ is dimensionless because 
it parameterizes the unit-size dimensionless function. We seek the best (where "best" will 
be defined below) ikf^^P- Gaussian MoG (where Af^^P is an integer) approximation 

Qexp(^) « 5]a;„^PiV(^|0,V-P) (6) 

m=l 

y-p ^ C'j , (7) 

where all of the means are exactly zero and all of the variances V^^p in the MoG can be 
represented as a scalar f^P multiplied by the identity matrix / because we are requiring 
this dimensionless function to be precisely circular (so every component is itself circular and 
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concentric). Similarly for the de Vaucouleurs (dev) profile 



a'^'''' = 7.66924944 (9) 
Q'i-(^) ^ 5^ ^^(^ I 0, (10) 

ydev dev J ^j) 



The half-light inverse-radius parameters a^^^ and a'^^^ are from Ciotti & Bertin (1999). The 
challenge we meet below is to determine the parameters 

l^m 5 Jm=l 5 l^m ) '^m Jm=l K^'^J 

to best approximate the traditional galaxy profile functions, under some sensible definition of 
the word "best" , as a function of the model complexity parameters (numbers of components) 
(M<'^P,M'i^^). 

In addition to these, there are general Sersic ("ser") profiles, of which the exp and dev 
profiles are special cases. The general ser profile has one parameter (the "index") n: 

gser(n)^^) _ exp(-a'^'=''(") [j^l^/" - 1]) (13) 

{a''''^^\a'^'^^\a'^''~^'>} = {3.67206075,5.67016119,9.66871461} , (14) 



where we have given the constant a^^'^^''^^ for just a few values of n ( Ciotti Sz Bertin||l999 the 
exp and dev profiles given above provide values for n = 1 and n = 4.) 



The SDSS pipelines ( Lupton eta/.||2001 ) make use of modified profiles, which have been 



truncated smoothly at large radius and (in the case of the de Vaucouleurs profile) "softened" 
at the center. The SDSS form of the exponential (lux) profile is 

f exp(-ai""[|^| - 1]) for 1^1 < 3 

Q'^^^(0 = <^ exp(-ai-[|^|-l]) [l-[|^|-3]f for 3 < 1^1 < 4 (15) 
[ for 4 < 1^1 

a'^^ = 1.67835 , (16) 
and the SDSS form of the de Vaucouleurs (luv) profile is 

f exp(-ai"^ + 0.0004]V8 _ ^j) |^| < 7 

Q'"^(|) = i exp(-ai"^ + 0.0004] V8- 1]) [1 - [1^1 -7]2]' for 7 < |^| < 8 (17) 
[0 for 8 < 1^1 

a^"^ = 7.66925 . (18) 
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The half-light inverse-radius parameters a^^^ and a^^"^ — and the softening and cutoff radius 
parameters — are taken from the SDSS codebase. 

The profiles above are normalized to have unit intensity (approximately) at their half- 
light radii. In many cases, the investigator wants profiles that are normalized to have unit 
total flux (intensity integrated over sohd angle) . Although there is an analytic result for the 

dev profile, numerical integration of the concentrated profiles dev and luv to determine total 
fluxes can be challenging. This is not true for the MoG approximations: Each Gaussian is 
normalized, so the sum of the amplitudes aj^^ (for the luv profile, say) gives the total 
flux for the MoG approximation to that profile. 

We seek the best MoG approximations. This necessitates definition of the word "best" . 
If we think of the profiles as being two-dimensional probability distribution functions (for, say, 
the arrivals of photons), then one natural choice is the K-L divergence or similar cross-entropy 
or information-theoretic measure. However, in typical astronomical imaging, the galaxy is 
superimposed on a substantial, flat sky level, and the noise in the data is close to Gaussian. 
This suggests more chi-squared-like objectives. We adopt the latter, in part because they 
are most appropriate for our speciflc proposed application (modeling SDSS-like astronomical 
imaging), but experiments we have performed suggest that information-theoretic objectives 
also lead to good results. 

In detail, the chi-squared objective we minimize — the badness — is a squared residual 
between the exact proflle function Q{^) and its MoG approximation. It is designed to 
be equivalent to a chi-squared statistic in a homoskedastic two-dimensional image of the 
proflle taken with extremely high angular resolution (pixels of size 0.001 the half-light radius) 
and vanishing point-spread function. Quantitatively the badness is defined to be the mean 
squared residual in the Q{$,) functions, which are normalized to have unit intensity at the 
half-light radius, averaged over a two-dimensional circular region in the ^ plane centered 
on the (circularly symmetric) profile and extending out to radius .^max- We use .^max = 8 
for all profiles except the lux profile, for which we use ^max = 4. In practice, the badness 
is computed in a one-dimensional numerical integral but the integral is weighted in radius 
(weight increasing linearly with radius) to make it equivalent to the two-dimensional chi- 
squared. We also add to the badness a very tiny coefficient (on the order of 10^'^ of the 
best-fit badness) times the sum of the variances Vm for regularization. In practice, this term 
doesn't have much effect and could be dropped. 

Optimization (minimization) of the badness is performed by the scipy implementation 
of the BFGS algorithm, with many initializations to explore multiple local minima. Further 
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details are available in the code, which is publicly availablejj 

The results of the optimizations are shown in Tables [T] and [2] and Figures [T] through [5] 
All the results shown in these figures and tables and more are available in machine-readable 
form from author DWH upon request. In the Tables and Figures we show root- variances 
rather than variances because these have units of half-light radii; they are simple standard 
deviations for the Gaussian components. 

In our work on The Tractor, we use the M^^^ = 6 lux and the M^"^ = 8 luv profiles. 
Our advice to users would be to do the same. We use the lux and luv over the exp and dev 
partly because of their better behaviors numerically, and partly because they and we are 
both part of the SDSS tradition. These — M^"^ = 6 and M'*^^ = 8 — are good compromises 
between mixture complexity (M) and quality of fit (badness). Also, even the best-fitting 
late-type and early-type galaxies deviate from exponential and de Vaucouleurs fits by more 
than do these high-quality MoG approximations; no precision is lost. 

In Figure [ij we show the dependence of amplitudes am^^""^ and variances Vm^^"''^ on the ser 
index n. There is clearly continuity; a valuable follow-up project would be to give expressions 
for the amplitudes and variances as a function of ser index n. In the absence of cleverness 
our advice would be to make use of smooth interpolation. 

The value of these MoG approximations comes when they are to be convolved with a 
PSF (usually in fact a pixel-convolved PSF) that is itself also represented as a MoG. In this 
scenario, the PSF iIj{Ax) — which is thought of as a function of focal-plane displacement Ax 
away from, say, a true stellar position — is represented as a i^-Gaussian MoG 

K 

ij{Ax) = J2pkN{Ax\mk,Vk) (19) 

k=l 
K 

1 = Y.P'^ ' (20) 

k=l 

where the means are not required to vanish because the PSF can have arbitrarily non- 
trivial structure (think speckles and the like) and the variances Vk will not in general be 
proportional to the identity or even diagonal because the PSF will not in general be round. 
An example that illustrates the use of this PSF is the following: A star of fiux Sg at focal- 
plane position Xg will lead to an image (PSF-convolved intensity map) of the form 

K 

I {x\ star, Ss,Xs) = ^ SsPk N{x\xs + mk,V k) ■ (21) 

k=l 



^ https : / / github . com/ davidwhogg/TheTractor/ 
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That is, when the PSF is represented as a MoG, any image of a star — or indeed any image 
of any set of stars — is also represented MoG. 

Applying this PSF to an exp or dev galaxy is slightly more complicated, because the 
galaxy has not just a flux Sg and a central position Xg] it also has a shape. Because 
we are only considering these simple galaxies, we are only permitting ellipsoidal shapes, 
which can be represented by a semi-major axis a, a semi-minor axis b, and a position angle 
(j), or equivalently by eigenvalues a,b and eigenvectors Ui,U2, or equivalently by an affine 
transformation Rg that takes a circle to the relevant ellipse (and is therefore a general 
representation of an ellipse; it is also the matrix square root of the symmetric variance 
tensor describing the ellipse). The galaxy is distorted by this affine transformation prior to 
PSF convolution, so the focal-plane image (PSF-convolved intensity field) for a general (say) 
exp galaxy is given by 

K M^^P 

I{X\ eXp,Sg,Xg,Rg) = ^ ^ S g tt^^^ P k N {x \ X g + Ulk , V gm + V k) (22) 

k=l m=l 

Vgm = Rg-VZ^-Rg" (23) 

Rg = [aui,bu2] , (24) 

where a and b are the major and minor axis lengths of the galaxy ellipse (in appropriate 
units) and Ui and U2 are the eigenvectors in image coordinates pointing in the major-axis 
and minor-axis directions respectively. Implicitly all vectors are two-dimensional column 
vectors, and Rg is a 2 x 2 affine transformation matrix that contains the "shape" (position 
angle, major-axis, and ellipticity) information about the galaxy. The dev, ser(n), luv, and 
lux cases are all essentially the same. Note the important and key result of this Note, to 
wit, that a MoG galaxy model (with M components) convolved with a MoG PSF model 
(with K components) yields a MoG model image (with [Mi^'] components). In Figure |6] we 
show how we are using these MoG approximations in The Tractor — a generative modeling 
framework for measuring astronomical objects — to render PSF-convolved galaxy images. 

In the above we said "pixel-convolved PSF". In every context, when modeling images, 
it is valuable to use the pixel-convolved PSF. With this definition of the PSF, the pixelized 
image is the PSF-convolved true model evaluated at the pixel centers. This operation is fast. 
Other definitions for the PSF (the non-pixel-convolved, for example) require that the user 
do two convolutions, the first with the PSF and the second with the square (or worse) pixel. 
Our advice: Only fit for and use pixel-convolved PSFs. 

If your PSF is not in MoG form, it is still the case that convolution of a MoG ap- 
proximation of a dev (say) profile will in general be easier than convolution of the original 
dev profile. The reason is that convolution of a Gaussian with any PSF is fast (indeed most 
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image-processing languages have such functions built in) ; the PSF-convolved profile becomes 
in this case just a mixture of Gaussian-convolved PSFs. 

The speed-ups that can be obtained by using MoG approximations can be very large. In 
our image-modeling project The Tractor, we were PSF-convolving by rendering the profiles 
(especially the profile centers) at very high resolution (hundreds to thousands of resolution 
elements in the central pixel are necessary for good precision on the dev profile). We were 
then convolving that high-resolution model with a low-resolution PSF and rendering to a 
low-resolution image pixel grid. These expensive operations were obviated by the MoG 
profiles, which involve only rendering a small number of Gaussians at the pixel centers on 
the low-resolution pixel grid. The MoG approximations saved us more than an order of 
magnitude in compute time, especially in optimization, where derivatives have to be taken 
with respect to the unconvolved model properties. 

In the SDSS, GALFIT, and much of our own work, the models that are fit are (ef- 
fectively) mixtures of exp and dev or exp and ser or lux and luv profiles. Mixtures of 
profiles that are each themselves mixtures of Gaussians are no harder to render than ei- 
ther profile separately. There is some book-keeping, of course, because each component gets 
affine-transformed separately before they are both PSF-convolved. 

One amusing aspect of MoG profiles has to do with projection from three to two dimen- 
sions. The projection of a three-dimensional Gaussian is a two-dimensional Gaussian; the 
two-dimensional, rigid, circular approximations we have made for the ser profiles can be de- 
projected to rigid, spherical approximations to the three-dimensional profiles trivially. The 
two-dimensional models we have started with are not accurate models of galaxies in detail — 
no galaxy follows exactly any ser profile — so deprojection of our approximations are not that 
interesting in themselves. However, the general program of fitting two-dimensional sources 
with MoGs may have strong implications in the future for three-dimensional modeling and 
deprojection. 

It is a pleasure to thank Brendon Brewer (Auckland), Jim Bosch (Princeton), and Kevin 
Bundy (IPMU) for valuable comments. This work was supported in part by NASA (grant 
NNX12AI50G) and the NSF (grant IIS-1124794). All the results shown in the figures and 
tables and more are available in machine-readable form from author DWH upon request. 
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Tabic 1: The amplitudes and root-variances for the best mixture-of-Gaussian approximations 
to the exp and dev profiles, for different mixture sizes. The total (dimensionless) fluxes and 
the badnesses are given for each approximation. 
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Table 2: Same as Table [T] but for the lux and luv profiles. 




Fig. 1. — top-left: The true exp profile (thin black line), the best JVP^^^ = 6 mixture- 
of-Gaussian approximation (thick grey line), and the component Gaussians (multiplied by 
their corresponding amplitudes) contributing to the approximation (thin grey lines). The 
plot title text gives ^max and the badness, top-right: The same but shown logarithmically. 
bottom-left: A representation of the residual or devation, on which the badness is computed. 
bottom-right: The same but shown fractionally and logarithmically. 
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dev / il/'^'^' = 10 / Cmax = 8 / badness = 3.90 x lO^^ dev / A/''"' = 10 / ^max = 8 / badness = 3.90 x lO"'^ 




dimensionless angular radius ^ dimensionless angular radius ^ 




Fig. 2. — The dev profile and the best M = 10 approximation. The panels are equivalent 
to those in Figure [TJ 
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lux / = 6 / ^max = 4 / badness = 4.64 x lO-** lux / il/''" = 6 / Cmax = 4 / badness = 4.64 x lO"*' 




dimensionless angular radius ^ dimensionless angular radius ^ 



luv / Ml"" = 8 / Uax = 8 / badness = 8.42 x IQ-^ luv / Af'" = 8 / Uax = 8 / badness = 8.42 x IQ-*' 




dimensionless angular radius ^ dimensionless angular radius ^ 



Fig. 3. — The lux and luv profiles and approximations. The panels are equivalent to those 
in the top-row of Figure [l] 
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ser2 / M""-- = 10 / Cmax = 8 / badness = 4.37 x 10"' 




dimensionless angular radius f 



ser3 / A/'"'-' = 10 / f„iax = 8 / badness = 8.38 x IQ-'' 




dimensionless angular radius ^ 



ser5 / M"*"^^ = 10 / gmax = 8 / badness = 1.02 x 10 approximations to ser profiles at AP'"' = 10 




dimensionless angular radius ^ root- variance ^/rf^'" (units of half-light radii) 



Fig. 4. — Three ser profiles — with n = 2, 3, and 5 — and approximations. The top-left, top- 
right, and bottom-left panels are equivalent to those in the top-right of Figure [T} bottom- 
right: The dependence on the amplitudes am^^^^ and root-variances "J^u^^^ on ser index 
n. 
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approximations to exp approximations to dev 




approximations to lux approximations to luv 




Fig. 5. — Comparisons of approximations, top-left: The dependence of amplitude and 
root-variance Vvm^ on M^^p for the exp profile, top-right: The same but for the dev profile. 
bottom-left: The same but for the lux profile, bottom-right: The same but for the luv 
profile. 





Fig. 6. — Demonstration of use of the profiles, or the imphcit generative model in this 
Note, top-left: The circular dimensionless M^"^ = 8 mixture-of-Gaussian approximation to 
the luv profile, represented on a very fine pixel grid, top-right: The ellipse representing 
the non-trivial affine transformation to be applied to the circular, dimensionless profile. 
middle-left: The sheared profile, middle-riglit: A K — 3 mixture-of-Gaussian model of the 
pixel-convolved point-spread function, represented on the very fine pixel grid, bottom-left: 
The sheared profile convolved with the PSF, represented on the very fine pixel grid, bottom- 
right: The sheared luv convolved with the PSF, but now shown on a reahstic pixel grid. 
Because by assumption the PSF is a pixel-convolved PSF, the representation on the realistic 
grid is found simply by interpolating to the pixel centers the mixture-component Gaussians. 



